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Abstract 


A  new  rate-dependent  cohesive  zone  model  for  use  in  impact  and/or  hydrodynamic  ram 
delamination  problems  is  presented.  The  traction  based  model  includes  a  damage- 
dependent  term  and  two  cohesive  zone  viscosity  parameters.  The  first  viscosity 
parameter  adjusts  the  yield  traction  of  the  material,  while  the  second  parameter  adjusts 
the  peak  or  maximum  traction.  This  new  rate-dependent  cohesive  model  is  implemented 
in  LS-DYNA  (971  beta  revision  5397),  an  explicit  time  integration  FEA  code,  by 
defining  a  user  cohesive  material  model.  It  may  be  used  with  shell  or  brick  elements. 
Previously  published  Double  cantilever  beam  (DCB)  experiments  of  epoxy  and 
PEEK/carbon-fiber  composites  are  modeled  in  order  to  validate  the  rate-dependent 
cohesive  model.  Also,  in  order  to  provide  a  methodology  to  obtain  data  to  determine 
material  constants  for  this  rate-dependent  cohesive  zone  model,  a  simple  inexpensive 
experimental  apparatus  for  high  rate  DCB  (Double  Cantilever  Beam)  sample  testing  was 
developed. 
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Preface 


A  new  rate-dependent  model  for  fracture  of  composite  materials  is  presented. 
Additionally,  an  inexpensive  experimental  apparatus  is  presented  which  provides  the  data 
necessary  for  the  rate-dependent  model.  Code  is  presented  which  can  be  used  as  a  user 
material  with  the  finite  element  explicit  package  LS-DYNA  (971  beta  revision  5397). 
This  new  rate-dependent  model  can  be  used  with  any  application  where  rate-dependent 
failure  is  important.  Specifically,  the  model  is  applied  to  survivability  of  composite 
material  under  hydrodynamic  ram.  The  model  presented  includes  only  mode  I  rate- 
dependence  and  may  be  used  with  shell  and  solid  elements. 

The  authors  thank  the  Aerospace  Survivability  and  Safety  Flight  of  the  46lh  Test  Wing  at 
Wright  Patterson  Air  Force  Base  for  continued  support  and  use  of  test  facilities. 


1  Summary 


A  new  rate-dependent  cohesive  zone  model  for  use  in  impact  mode  I  delamination 
problems  is  presented.  The  traction  based  model  includes  a  damage-dependent  term  and 
two  cohesive  zone  viscosity  parameters.  The  first  viscosity  parameter  adjusts  the  yield 
traction  of  the  material,  while  the  second  parameter  adjusts  the  peak  or  maximum 
traction.  LS-DYNA  (971  beta  revision  5397),  an  explicit  time  integration  FEA  code,  is 
used  to  model  double  cantilever  beam  (DCB)  experiments  of  fiber  composites  with 
loading  rates  varied  from  0.65  m/s  to  20  m/s.  The  rate-dependent  cohesive  zone  model  is 
used  to  formulate  a  rate-dependent  cohesive  element  for  use  with  shell  or  solid  elements. 
This  rate-dependent  cohesive  model  is  implemented  in  LS-DYNA  by  defining  a  user 
cohesive  material  model.  Using  previously  published  data,  epoxy  and  PEEK/carbon-fiber 
composites  are  found  to  require  a  rate-dependent  cohesive  model  to  adequately  model  the 
experiments.  The  rate-dependent  cohesive  zone  model  presented,  with  a  single  set  of 
constants  for  each  material,  produces  agreement  between  the  experimental  and  FEA 
results. 

A  simple  inexpensive  experimental  apparatus  for  high  rate  DCB  (Double  Cantilever 
Beam)  sample  testing  was  developed.  The  apparatus  uses  a  cocked  spring  mechanism  to 
provide  high  rate  mode  I  loading  to  a  pre-cracked  DCB  specimen.  Testing  of  six  samples 
was  performed.  Loading  rates  varying  from  0.6  to  14.5  m/s  were  obtained.  A  finite 
element  analysis  simulating  the  testing  of  DCB  samples  was  performed.  Comparisons 
with  the  rate-independent  (v=//=0)  and  rate-dependent  (rtf 0)  cohesive  model  and  the 
experimental  data  were  performed.  It  was  determined  that  the  rate-dependent  model 
provided  better  agreement  with  the  experimental  data  than  the  rate-independent  model. 
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Introduction 


2.1  Background 

The  use  of  polymer  matrix  composites  and  bonded  and/or  co-cured  assemblies  in 
airframe  structures  has  shown  promise  in  achieving  the  performance  and  cost  goals  of 
next  generation  fighter/attack  aircraft.  The  weight  and/or  affordability  benefits  may  be 
limited,  however,  by  the  need  to  meet  survivability  requirements.  The  current  survivable 
design  procedure  is  to  size  a  structure  for  flight,  fuel  pressure,  crash,  etc.  loads  and  then 
ballistically  test  the  resulting  design  to  determine  its  survivability  capability.  For  metal 
structures,  this  remains  a  feasible  process  since  there  are  plenty  of  historical  ballistic  test 
data  available  for  use  in  developing  design  requirements.  However,  this  is  not  the  case 
for  composite  structures. 

The  typical  manned  fighter  survivability  design  goal,  defined  by  the  Live  Fire  Law[l],  of 
55%  design  limit  load  (DLL)  residual  strength  after  damage  from  a  30mm  high  explosive 
incendiary  (HEI)  threat  has  not  yet  been  demonstrated  on  an  all-composite  platform. 
(Unmanned  vehicles  and  non-fighter  aircraft  will  have  less  stringent  requirements.)  In 
attempting  to  meet  the  Live  Fire  Law  requirements,  the  F-22  program  conducted  several 
ballistic  tests  on  all-composite  wing  designs  without  success.  Ballistic  tolerance  was 
finally  accomplished  by  replacing  five  bolted  composite  spars  with  five  bolted  titanium 
spars.  This  cost  the  program  thousands  of  dollars  per  aircraft  and  schedule  slippage  and 
development  costs  as  well  as  additional  weight  in  each  aircraft. 

The  Air  Force  and  Navy  consider  the  survivability  problem  of  an  all-composite  structure 
critical  and  have  invested  a  significant  amount  of  funding  to  address  this  issue.  The  Air 
Force’s  Decoupled  Fuel  Cells  (DFC)  program  [2]  identified  the  benefits  of  addressing 
survivability  earlier  in  the  design  phase  with  cost  and  weight  savings  using  a  cocured 
wing  design  for  the  F-22  and  eliminating  the  titanium  spars.  This  study  demonstrated 
that  a  wing  design  that  relies  on  a  bolted  metal  substructure  to  meet  the  live  fire 
requirements  for  combat  aircraft  can  cost  and  weigh  substantially  more  than  a  survivable 
composite  design.  Furthermore,  the  promise  of  future  weight  and  cost  savings  can  only 
be  realized  by  addressing  the  survivability  of  composite  structures  early  in  the  design 
phase  of  an  aircraft. 

Advanced  processing  techniques,  interlaminar  reinforcement  technologies  and  innovative 
design  concepts  have  been  developed  in  recent  years  and  provide  significant 
improvements  towards  achieving  a  survivable,  all-composite  wing  structure  while 
minimizing  its  weight  and  cost.  Until  the  execution  of  the  Composites  Affordability 
Initiative  (CAI)  in  the  late  1990’s,  industry  lacked  a  set  of  rules  and  procedures  that  could 
guide  engineers  in  incorporating  these  survivability  features  early  in  the  structural  design 
process.  The  Survivability  Group  of  the  CAI  Pervasive  Team  drafted  the  HRAM 
Survivability  Design  Guidelines  [3],  which  contained  this  direction  so  that  the  weight  and 
cost  goals  could  be  met.  The  use  of  these  guidelines  will  reduce  the  design  and  analysis 
development  time  by  greatly  minimizing  the  need  for  redesigns  after  the  ballistic  testing 


2 


is  completed.  Their  use  will  also  potentially  reduce  the  amount  of  costly  testing  needed 
by  providing  a  history  of  ballistic  test  results  and  proven  designs  for  damage 
containment.  Ballistic  testing  cannot  be  completely  eliminated  from  the  process,  but  the 
cost  of  testing  can  be  reduced  by  testing  smaller  articles  and  incorporating  historical  test 
data  into  the  design  process. 

The  general  sequence  for  designing  for  survivability  (as  specified  in  the  Design 
Guidelines)  is  as  follows: 

1)  Develop/Identify  the  survivability  requirements  for  the  program , 

2)  Design  the  structure  for  flight  and  hydrodynamic  ram  (HRAM)  loads,  conduct 
developmental  testing,  and 

3)  Ballistically  test  the  final  designs,  either  full-scale  articles  or  large 
subcomponents. 

Implied  within  2)  above  is  the  assumption  that  the  designer  has  validated  and  reliable 
modeling  and  simulation  (M&S)  tools  at  his  disposal  in  order  to  perform  analyses  which 
will  assist  him  in  converging  on  an  optimum  design.  Although  the  Survivability  Group 
of  the  CAI  Pervasive  Team  made  significant  progress  in  this  area,  they  identified  two 
shortcomings  specifically  with  the  hydrodynamic  ram  analysis  techniques,  which  must  be 
resolved  before  the  tools  can  be  declared  robust.  These  shortcomings  are: 

1)  More  work  is  needed  to  accurately  model  the  details  that  differentiate  one  joint 
design  from  another.  The  incorporation  of  stitching,  z-pinning,  co-curing,  and 
bonding  of  the  joints  leads  to  different  failure  modes  and  paths.  The  area  in 
which  most  work  is  needed  is  in  the  proper  selection  of  failure  criteria  within  the 
joints  as  well  as  the  entire  model.  Right  now,  the  state-of-the-art  is  to  use  the 
elastic-plastic  smeared  properties  technique  in  modeling  the  structure  and  “fuze 
elements”  to  model  the  joints. 

2)  High  loading  rates  encountered  in  hydrodynamic  ram  (HRAM)  events  are 
neither  well  understood  nor  well  characterized.  As  rate  dependent  material 
properties  are  determined  and  documented,  the  analysis  codes  must  have  the 
capability  to  incorporate  them.  Although  there  is  limited  capability  in  this  regard, 
more  work  is  needed. 

These  issues  are  significant  in  accurately  predicting  the  survivability  of  multi-sparred, 
fuel-filled  composite  wing  structures.  The  plan  for  such  wings  is  to  have  enough 
structural  integrity  remaining  after  impact  to  carry  the  loads.  If  the  design  tools  aren’t 
accurate,  they  may  over  or  under  predict  the  viability  of  the  remaining  structure.  In  either 
case,  increased  cost  (in  terms  of  dollars,  time,  and  human  life)  may  ultimately  result. 

Historically,  the  failure  of  joints  is  the  result  of  fast  fracture  that  occurs  in  resin  rich 
zones  of  the  composite  joints.  The  development  and  implementation  of  fracture 
mechanics  concepts  into  explicit  time  integration  FEA  coding  is  required  if  such  events 
are  to  be  modeled.  Additionally,  the  development  of  a  fracture  mechanics  based  model 
for  energy  release  rate  as  a  function  of  the  crack  opening  rate  is  needed  to  account  for  the 
high  loading  rates  caused  by  HRAM. 
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2.2  Project  Objectives 

There  were  three  main  objectives  of  the  proposed  work: 

1)  Develop  a  fracture  mechanics  based  model  for  energy  release  rate  as  a  function  of 
crack  opening  displacement  (COD)  rate  in  co-cured  and  bonded  organic 
composites. 

a)  Conduct  basic  research  into  improving  the  current  stress/strain  failure  criteria 
which  are  not  robust  enough  to  model  high  loading  rate  fracture  failure. 

b)  Develop  a  basic  method  for  acquiring  energy  release  rate  data  that  is  high 
loading  rate  dependent.  A  simple  double  cantilever  beam  (DCB)  specimen  will  be 
used  for  the  data  collection. 

c)  Use  the  data  from  l.b  in  order  to  develop  an  empirical  based  model  relating 
COD  rate  to  energy  release  rate.  Using  fracture  mechanics  concepts  a  release 
pressure  (traction-separation  law)  will  then  be  derived. 

2)  Develop  a  finite  element  formulation  for  incorporating  into  an  explicit  time 
integration  finite  element  analysis  package  a  nodal/element  based  method  for 
describing  a  force/stress  based  fracture. 

a)  Incorporate  the  best  practices  found  in  their  basic  research  of  la.  through  lc. 
into  a  finite  element  formulation. 

b)  Demonstrate  through  modeling  and  simulation  of  a  simple  DCB  specimen  the 
applicability  of  the  new  high  loading  rate  failure  mechanism. 

3)  Validate  the  new  finite  element  formulation  and  high  COD  rate  fracture 
mechanics  based  failure  mechanism. 

a)  RHAMM  Technologies,  LLC  will  complete  validation  by  comparing  the  DCB 
modeling  results  with  DCB  experimental  high  strain  rate  data. 

b)  RHAMM  Technologies,  LLC  will  build  a  finite  element  model  of  simple 
composite  wing  box  and  demonstrate  joint  failure  using  the  new  failure 
mechanism  under  HRAM  loading. 

All  objectives  with  exception  of  3b  were  met.  The  personnel  involved  in  this  work  are 
presented  in  Appendix  Al.  The  publications  resulting  from  this  work  are  presented  in 
Appendix  A2. 

2.3  Statement  of  scope 

All  Live  Fire  Testing  and  Evaluation  (LFT&E)  and  Joint  Live  Fire  (JLF)  and  other 
vulnerability  assessments  are  negatively  impacted  by  the  lack  of  analysis  capability. 

This  is  particularly  crucial  to  the  design  of  new  aircraft,  such  as  the  Joint  Strike  Fighter 
(JSF),  which  have  multiple  live  fire  requirements.  Because  the  JSF  and  future  aircraft 
will  use  significantly  more  composites  than  their  predecessors,  new  design  tools  are 
needed.  Without  a  complete  analysis  capability  that  includes  the  ability  to  model 


4 


composite  joint  concepts,  designers  are  forced  to  build  and  test  composite  boxes  that  are 
extremely  expensive.  When  fully  developed,  this  concept  will  provide  tools  for 
accurately  predicting  the  survivability  of  multi-sparred,  fuel-filled  composite  wing 
structures. 

2.4  Report  Organization 

The  experimental  apparatus  for  testing  DCB  samples  at  high  strain  rates  is  presented  in 
Section  3.1.  The  new  rate-dependent  cohesive  model  is  presented  in  Section  3.2,  while 
the  details  of  the  LS-DYNA  FEA  model  are  presented  in  Section  3.3.  The  validation  of 
the  rate-dependent  cohesive  element  developed  in  this  research  is  presented  in  Section 
4. 1 .  This  validation  compares  experimental  data  at  various  rates  of  composite  DCB 
samples  from  literature  to  the  FEA  results  using  cohesive  model  with  and  without  rate- 
dependence.  Section  4.2  presents  the  experimental  data  obtained  using  the  high  rate  DCB 
testing  apparatus  developed  during  this  study.  Finally  in  Section  4.3,  the  cohesive  model 
is  used  to  model  these  experiments  and  the  results  are  compared  to  the  experimental  data. 
Conclusions  are  presented  in  Section  5. 
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3  Methods,  Assumptions  and  Procedures 

3.1  High  Strain  Rate  Experimental  Apparatus 

An  experimental  apparatus  that  may  be  used  to  test  DCB  test  articles  at  high 
loading/strain  rates  was  designed,  built  and  tested.  This  section  describes  the  DCB 
specimens  testing,  the  experimental  apparatus  design,  and  the  experimental  results. 

3.1.1  DCB  Experimental  Specimen 

A  schematic  drawing  and  photograph  of  the  DCB  experimental  specimen  are  shown  in 
Fig.  1.  The  specimens  are  25.4  mm  wide,  229  mm  long,  and  3.05  mm  thick,  with  an 
initial  crack  of  25.4  mm  at  one  end,  artificially  simulated  by  a  Teflon  separator  25.4  mm 
wide  and  0.025  mm  thick,  inserted  in  the  mid-surface  of  the  laminate.  T-shaped 
aluminum  tabs  (12.7  mm  wide,  25.4  mm  long,  and  3.18  mm  thick)  are  attached  to  the 
DCB  specimens  using  a  two  part  epoxy  adhesive.  White  reference  dots  were  painted  on 
the  side  of  each  DCB  specimen  to  facilitate  crack  velocity  and  deflection  measurements. 


Applied  Load 


Deflection 


Figure  1.  DCB  experimental  specimen  setup  and  loading:  schematic  and  photograph. 
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3.1.2  DCB  Experimental  Apparatus 

The  experimental  testing  setup  includes  a  high  rate  loading  test  machine,  high  speed 
digital  video  camera  (Phantom  V7),  TEMA  (TrackEye  Motion  Analysis)  software,  and 
LabView  data  acquisition  system  for  collecting  strain  and  load  data.  The  drawings  of  this 
apparatus  are  contained  in  Appendix  A3. 

The  high  rate  loading  test  machine  was  designed  to  store  energy  in  a  compression  spring 
and  quickly  release  the  energy  to  load  a  double  cantilever  beam  specimen.  A  cable  is 
attached  to  each  tab  on  the  DCB.  Each  cable  provides  the  loading  as  shown  in  Fig.  2. 

One  end  of  the  cable  is  fixed  (load  cell  end)  and  one  end  is  connected  to  a  cable  with  an 
aluminum  ball  attached  to  the  other  end.  The  aluminum  ball  rests  on  a  compression 
spring  which  is  compressed  and  held  with  a  lever  until  a  trigger  is  pulled.  When  the 
trigger  is  pulled,  the  spring  accelerates  the  aluminum  ball  pulling  on  the  cable  and  pulling 
apart  the  DCB.  Fig.  2  shows  a  sketch  of  the  configuration  at  time  t=0  and  some  time  t>0 
after  the  trigger  has  been  pulled.  Cocking  the  spring  was  accomplished  by  using  a  hand 
actuated  hydraulic  actuator  to  pull  up  on  the  cable  attached  to  the  compression  spring 
before  the  specimen  is  loaded  into  the  apparatus. 


Figure  2.  DCB  test  machine  at  t=0  and  t>0 

3.1.3  Data  Collection  and  Analysis 

A  piezoelectric  dynamic  load  cell  is  used  to  determine  the  applied  load  on  the  fixed  end. 
Deflection  at  the  applied  load  and  crack  velocity  is  measured  from  high  speed  digital 
video  using  TEMA  software  and  white  reference  dots  on  the  side  of  the  DCB  specimen. 
The  Phantom  V7  digital  camera  is  set  to  collect  4800  frames  per  second  at  800x600 
pixels  resolution.  The  digital  video  along  with  the  TEMA  software  package  provides  the 
necessary  capability  to  extract  both  crack  tip  velocity  and  crack  opening  displacement 
(COD). 
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3.2  Rate  Dependent  Cohesive  Model  Formulation 

Cohesive  models  have  been  widely  used  to  model  dynamic  crack  growth  (Chowdhury 
and  R.  Narasimhan  2000;  Rahul-Kumar  et  al.  2000;  Ruiz  et  al.  2001;  La  Saponara  et  al. 
2002).  Many  models  have  been  introduced  including:  extrinsic,  intrinsic,  perfectly  plastic, 
linear  softening,  progressive  softening,  and  regressive  softening  (Camanho  and  Davila 
2004).  Several  rate-dependent  models  have  also  been  introduced  (Glennie  1971;  Xu  et  al. 
1991;  Costanzo  and  Walton  1997;  Tvergaard  and  Hutchinson  1996;  Fager  and  Bassani 
1991;  Kubair  and  Geubelle  2003).  A  rate-dependent  cohesive  zone  model  was  first 
introduced  by  Glennie  (1971)r  =  r/(l  +  tju) ,  where  the  traction  in  the  cohesive  zone  is  a 
function  of  the  crack  opening  displacement  time  derivative.  Xu  et  al.  (1991)  extended  this 
model  by  adding  a  linearly  decaying  damage  law  r  =  rc(l  -  u/ucr)(l  +  tj  u/cs ) .  In  each 

model  the  viscosity  parameter  (t|)  is  used  to  vary  the  degree  of  rate  dependence.  Kubair 
et  al.  (2003)  thoroughly  summarized  the  evolution  of  these  rate-dependant  models  and 
provided  the  solution  to  the  mode  III  steady-state  crack  growth  problem  as  well  as 
spontaneous  propagation  conditions. 

The  rate-dependent  cohesive  zone  model  resulting  from  this  research  is  presented  in  Eq. 
(1-2).  This  model  is  unique  it  contains  a  rate-dependent  term  for  the  onset  of  the  opening 
(wi )  of  the  cohesive  zone  combined  with  a  maximum  traction  rate-dependent  term  (utj  ). 
The  model  is  also  unique  in  that  a  cohesive  zone  starting  or  yielding  displacement  ( us ) 
may  be  defined  to  greatly  increase  numerical  stability.The  critical  strain  energy  release 
rate,  found  by  integrating  over  the  appropriate  crack  opening  displacements,  is  given  by 
Eq.  (3). 


r  = 


(r  +vu)— 


for u<u  , 


+W77-(r  +  VWt)  /  x 
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For  mode  II  and  III  loading  the  cohesive  model  is  rate-independent  as  presented  in 
equations  (4-6).  The  mode  II  and  III  crack  opening  displacements  are  combined 
according  to  equation  (4)  in  order  to  form  the  damage  constant  y/2i .  The  cohesive 
tractions  are  assigned  with  equations  (5  and  6).  This  is  a  simple  bi-linear  model.  The 
constant  K2i  is  used  for  both  mode  II  and  III  loadings  for  the  hardening  portion  of  the 
traction-separation  curve,  however  two  separate  constants  could  easily  be  defined  if 
needed.  A  rate-dependent  model  for  mode  II  and  III  is  left  to  another  research  work. 
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An  example  traction  separation  curve  for  the  mode  I  rate-dependent  portion  of  the  model 
is  shown  in  Figure  3.  The  initial  portion  of  the  traction  separation  curve  is  defined  by  two 
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parameters:  the  yield  stress  x y  and  the  displacement us .  Typically,  when  the  traction 
exceeds  x  ,  this  signals  the  nucleation  of  the  cohesive  zone.  However,  the  author’s  found 
it  necessary  for  us  to  be  greater  than  zero  but  small  (in  this  case,  0.001  pm)  in  order  to 

facilitate  numerical  stability  in  the  time  explicit  finite  element  numerical  methodology.  In 
the  case  of  finding  an  analytical  solution  of  a  problem  employing  this  cohesive  model, 
the  value  of  us  may  be  set  to  zero. 

The  traction  rmax  is  the  maximum  traction  that  can  occur  if  no  rate  dependence  (v  =77= 0) 

is  present.  This  traction  occurs  at  the  critical  displacement  ucr ■  For  crack  opening  values 
larger  than  ucr,  the  traction  reduces  linearly  to  zero  at  the  crack  opening  of  uend-  This 
occurs  through  the  damage  parameter  \j/  varies  from  one  (at  u=  ucr)  to  zero  (at  u-ue„d). 
Through  the  viscosity  term  v,  rate-dependence  is  given  to  the  yield  traction xy,  while  the 

viscosity  term  p  gives  rate  dependence  to  rmax .  The  traction  separation  curve  remains 
linear  between  us  and  ucr,  and  between  ucr  and  ue„d  regardless  of  the  values  of  v  and  77. 
Figure  4  presents  an  example  of  a  traction  separation  curve  for  an  arbitrarily  chosen 
displacement  function  u(t)=500t3.  For  this  example  rj  is  fixed  (77=0)  and  v  is  varied.  As 
the  cohesive  zone  viscosity  increases,  the  critical  strain  energy  release  rate  ( Gc ), 
represented  by  the  area  of  the  traction  separation  curve,  increases  linearly  (see  Figure  5). 
Figures  6  and  7  present  the  effect  of  varying  the  viscosity  parameter  77  (with  v  =0)  on  the 
traction  and  critical  strain  energy  release  rate,  respectively.  Figure  8  shows  the  effect  of 
varying  both  v  and  77  with  77  =2  v.  It  should  be  noted  that  negative  values  of  v  and  77  are 
permissible  with  the  criteria  thatr  >  0 .  Thus,  when  77  is  less  zero  and  u  is  great  enough  to 
cause  xmax  to  decrease  to  zero  or  below,  the  traction  would  necessarily  be  set  to  zero. 

This  model  implemented  in  LS-DYNA  as  a  cohesive  user  defined  material  property.  The 
code  contained  in  the  file  dyna2 1  .f  is  presented  in  Appendix  A4.  The  cohesive  element  in 
LS-DYNA  is  currently  only  available  in  a  beta  version  of  LS-DYNA  (971  beta  revision 
5397) 
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Figure  5.  Critical  strain  energy  release  rate  versus  v,  example  ior  armiraniy  cnosen  u(t)  wnn  vanea  v  (wun 
V=0)- 


Figure  6.  Cohesive  traction  separation  curves,  example  for  arbitrarily  chosen  u(t)  with  varied  rj  (with  v  =0). 
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3.3  LS-DYNA  FEA  Modeling 

The  DCB  specimen  geometery  modeled  in  this  study  is  presented  in  Figure  9.  The  DCB 
specimen  was  modeled  with  12500  shell  elements,  13052  nodes,  and  5500  cohesive 
elements.  The  specimen  or  shell  element  thickness  {hi  and  h2)  of  each  arm  of  the  DCB 
specimen  was  set  to  1.5  mm.  An  initial  crack  length  of  approximately  30  mm  was  used, 
depending  on  the  values  evident  from  the  reported  data.  The  reported  velocity  was 
prescribed  on  the  upper  arm  at  a  location  5  mm  from  the  pre-cracked  end  of  the 
specimen.  The  displacement  of  the  bottom  arm  was  fixed,  while  allowing  free  rotation  at 
the  same  location.  The  material  properties  obtained  by  Blackman  et  al.  (1995)  were  used. 
For  the  PEEK  composite  material,  the  bending  modulus,  density,  and  Poisson’s  ratio 
were  120  GPa,  1540  Kg/m3,  and  0.28,  respectively.  For  the  epoxy  composite  material,  the 
bending  modulus,  density,  and  Poisson’s  ratio  were  1 15  GPa,  1566  Kg/m3,  and  0.27, 
respectively.  The  effective  bending  modulus  was  found  using  a  3-point  bend  test 
(Blackman  et  al.  1995).  The  bending  modulus  was  found  to  be  rate-independent  using  an 
ultrasonic  test  rig  utilizing  Lamb  waves.  The  critical  portions  of  the  LS-DYNA  input 
(.key)  files  are  included  in  appendix  A5. 

The  crack  length  history  data  for  each  analysis  are  obtained  by  first  running  LS-Post.  In 
LS-Post  a  line  of  nodes  running  in  the  direction  of  crack  growth  are  selected  on  each 
cross  section  (upper  and  lower).  The  displacements  in  the  direction  of  interest  are  then 
saved  to  a  text  file.  The  crack  history  extraction  program  given  appendix  A6  is  then  run. 
The  user  is  prompted  for  filename  where  the  displacement  data  is  contained,  the  initial 
crack  length,  and  the  filename  that  contains  the  node  XYZ  coordinates. 
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Figure  9.  Finite  element  model  of  Blackman  et  al.  (1995;  1996)  DCB  sample. 
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Results  and  Discussion 


4. 1  Rate  Dependent  Cohesive  Model  Results  Validation 

The  rate-dependent  cohesive  zone  model  (section  3.2)  was  implemented  using  a  user 
defined  cohesive  material  model  in  LS-DYNA,  an  explicit  time  integration  FEA  code. 
The  experiments  performed  by  Blackman  et  al.  (1995;  1996)  were  modeled.  Two 
materials  were  analyzed  PEEK/carbon-fiber  and  Epoxy/carbon-fiber  composites  at  the 
testing  rates  reported  by  Blackman  et  al.  (1995;  1996)  varing  from  0.65  m/s  to  20.5  m/s. 
The  specimen  geometry  and  boundary  conditions  used  in  the  FEA  model  are  presented  in 
Section  3.3,  while  the  rate-independent  (v=  rj-0)  and  rate-dependent  FEA  results  are 
presented  in  Sections  4.3.1  and  4.3.2,  respectively. 

4.1.1  FEA  Rate-Independent  Results 

The  crack  length  (a)  versus  time  curves  reported  by  Blackman  et  al.  were  compared  to 
the  rate-independent  (v=>/=0)  cohesive  zone  FEA  results  for  the  appropriate  loading  rates. 
The  critical  strain  energy  release  rate  (Gc),  according  to  Eq.  (3),  was  set  so  that  agreement 
between  the  experimental  and  finite  element  analysis  crack  length  (a)  versus  time  curves 
was  obtained.  Example  crack  length  versus  time  plots  for  the  PEEK  and  epoxy  composite 
materials  are  presented  in  Figs.  10-13;  showing  a  comparison  of  the  rate-independent 
FEA  results  and  the  experimental  results.  Excellent  agreement  was  found  between  the 
rate-independent  and  experimental  results;  however  adjustment  of  the  Gc  was  required  for 
each  rate  in  order  to  reach  this  agreement.  Blackman  et  al.  observed  unstable  “stick-slip” 
crack  growth  of  the  PEEK  composite  at  a  rate  of  1.1  m/s.  The  rate-independent  model 
(v=ri=Q)  was  not  able  to  capture  this  behavior.  Gc  tended  to  increase  with  increasing 
loading  rate  for  the  epoxy  composite  material  as  shown  in  Fig.  14,  however  for  the  peek 
composite  material  Gc  decreased  with  loading  rate.  In  Fig.  14,  the  values  of  Gc  obtained 
for  this  study  are  compared  with  those  reported  by  Blackman  et  al.  (1995).  Reasonable 
agreement  was  found  considering  lack  of  complete  information  regarding  the  testing 
methodology.  Table  1  presents  the  tabulated  Gc  values  required  for  each  case.  From  these 
results  it  may  be  concluded  that  a  rate-dependent  cohesive  zone  model  is  required  to 
adequately  model  these  materials  and  experimental  results  when  significant  variation  in 
loading  rate  is  expected. 
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Figure  10.  PEEK/carbon-fiber  composite  FEA  and  experimental  (Blackman  et  al.  1995;  1996)  crack  length 
versus  time  for  a  loading  rate  1.1  m/s. 


Figure  11.  PEEK/carbon-fiber  composite  FEA  and  experimental  (Blackman  et  al.  1995;  1996)  crack  length 
versus  time  for  loading  rates  of  6.5, 14.9,  and  18.4  m/s. 
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Figure  12.  Epoxy/carbon-fiber  composite  FEA  and  experimental  (Blackman  et  al.  1995;  1996)  crack  length 
versus  time  for  loading  rate  of  0.65  m/s. 


Figure  13.  Epoxy/carbon-fiber  composite  FEA  and  experimental  (Blackman  et  al.  1995;  1996)  crack  length 
versus  time  for  loading  rates  of  8.0, 15.0,  and  20.5  m/s. 
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Loading  Rate,  m/s 


Figure  14.  Critical  strain  energy  release  rate  (Eq.  (3))  versus  loading  rate  of  rate-independent  model  (v 
=>7=0)  compared  to  reported  experimental  values  (Blackman  et  al.  1995;  1996). 


Table  1.  Tabulated  values  of  critical  strain  energy  release  rate  as  a  function  of  loading 
rate  for  the  current  rate-independent  model  (v  =rj  =  0). 


Epoxy/Carbon-Fiber 

Peek/Carbon-Fiber 

Loading  Rate , 
m/s 

Strain  Energy  Release 
Rate,  N/mm 

Loading  Rate , 
m/s 

Strain  Energy 
Release  Rate,  N/mm 

0.65 

0.168 

1.1 

1.241 

8.0 

0.251 

6.5 

0.897 

15.0 

0.141 

14.9 

0.732 

20.5 

0.278 

18.4 

1.103 

4. 1 .2  FEA  Rate-Dependent  Results 

An  FEA  was  performed  using  the  rate-dependent  cohesive  zone  model  (v£0  or  rpjtO).  The 
viscosity  terms  v  and  tj  and  the  crack  opening  parameters  ucr  and  ue„d  were  adjusted  until 
agreement  between  the  FEA  and  experimental  crack  length  versus  time  curves  for  the 
rates  analyzed  in  this  study  was  found.  A  single  set  of  parameters  was  desired  for  all 
loading  rates  analyzed.  Table  2  presents  the  rate-dependent  cohesive  element  parameters 
obtained  for  the  epoxy  and  PEEK  composite  materials.  As  the  trend  for  the  PEEK 
composite  material  was  a  decreasing  Gc  with  increasing  rate,  a  negative  value  of  tj  was 
required.  The  viscosity  parameter  v  was  chosen  to  be  zero,  as  in  this  particular  DCB 
configuration  adjusting  its  value  did  not  have  a  significant  effect  on  the  overall  crack 


18 


propagation  behavior.  Figs.  15  and  16  present  the  crack  length  versus  time  curves  at  each 
of  the  loading  rates  analyzed  for  the  epoxy  composite.  Good  agreement  was  found 
between  the  rate-dependent  FEA  and  the  experimental  results  for  loading  rates  ranging 
between  0.65and  20.5  m/s.  Figures  17  and  1 8  present  the  results  for  the  PEEK  composite 
for  rates  varied  from  1.1  to  18.4  m/s.  Good  agreement  between  the  FEA  and 
experimental  results  was  found.  Interestingly,  the  “stick-slip”  crack  propagation  behavior 
of  the  PEEK  composite  material  at  a  loading  rate  of  1 . 1  m/s  was  captured  with  this  rate- 
dependent  cohesive  model.  The  behavior  is  also  observed  at  the  higher  rate  of  18.4  m/s  as 
also  evidenced  in  the  experimental  results.  Thus,  with  the  rate-dependent  cohesive  zone 
model,  a  single  set  of  material  constants  for  each  material  was  able  to  model  the  rate- 
dependent  behavior  of  the  epoxy  and  PEEK  composites. 


Table  2.  Rate-dependent  model  cohesive  element  parameters. 


Parameters 

Epoxy/Carbon-Fiber 

PEEK/Carbon-Fiber 

uSt  Jim 

0.01 

0.01 

Mcr.M m 

0.344 

70.0 

M  pm 

7.0 

100.0 

Ty ,  GPa 

0.0200 

0.0200 

r  ,  GPa 

0.0275 

0.0275 

v,  pm/ms 

0.0 

0.0 

rj,  pm/ms 

50.0 

-5.30 

Figure  15.  Rate-dependent  model  results:  epoxy/carbon-fiber  composite  FEA  and  experimental  (Blackman 
et  al.  1995;  1996)  crack  length  versus  time  for  loading  rate  of  0.65  m/s. 
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Figure  16.  Rate-dependent  model  results:  epoxy/carbon-fiber  composite  FEA  and  experimental  (Blackman 
et  al.  1 995;  1 996)  crack  length  versus  rime  for  loading  rates  of  8.0,  1 5.0,  and  20.5  m/s. 


Time,  ms 

Figure  17.  Rate-dependent  model  results:  PEEK/carbon-fiber  composite  FEA  and  experimental  (Blackman 
et  al.  1995;  1996)  crack  length  versus  time  for  loading  rate  of  1.1  m/s. 
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Figure  18.  Rate-dependent  model  results:  PEEK/carbon-fiber  composite  FEA  and  experimental  (Blackman 
et  al.  1995;  1996)  crack  length  versus  time  for  loading  rates  of  6.5, 14.9,  and  18.4  m/s. 


4.2  Experimental  Results  from  High  Strain  Rate  Apparatus 

The  experimental  apparatus  as  described  in  section  3.1.2  was  used  to  test  the  DCB 
experimental  specimens  described  in  3.1.1.  The  DCB  specimens  composed  ofNCT-350- 
TR-50,  obtained  from  Newport  Adhesives  &  Composites,  Inc. 

The  crack  opening  displacements  along  the  direction  of  crack  growth  and  load  were 
collected  for  6  specimens.  Representative  results,  crack  length,  crack  velocity,  and  load 
time  histories,  are  shown  in  Figures  19-24.  The  testing  procedure  was  found  to  produce 
loading  rates  ranging  from  0.6  to  14.5  m/s. 
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Figure  19.  DCB  experimental  results  samples  1, 2, 4:  crack  length  versus  time. 
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Figure  21.  DCB  experimental  results  for  samples  1,  2,  and  4:  crack  velocity  versus  time. 


Figure  22.  DCB  experimental  results  for  samples  5,  6,  and  9:  crack  velocity  versus  time. 
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Figure  23.  DCB  experimental  results  for  samples  2  and  4:  measured  load  versus  time. 


Time,  ms 


Figure  24.  DCB  experimental  results  for  samples  5,  6,  and  9:  measured  load  versus  time. 
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4.3  Finite  Element  Analysis  of  Experimental  Results 

A  finite  element  analysis  simulating  testing  of  the  6  DCB  samples  (Samples  1,  2, 4,  5,  6, 
and  9)  was  performed.  The  measured  crack  mouth  opening  displacements  (CMOD)  as  a 
function  of  time  are  presented  in  Figures  25-26.  Piecewise  linear  curve  fits  were 
performed,  in  order  to  obtain  the  prescribed  velocities  for  the  model  loading  condition. 
The  velocities  used  are  shown  in  Table  3.  The  rate-independent  cohesive  model  (v=/p 0) 
was  used  with  a  single  set  of  parameters  presented  in  Table  4.  These  results  are  presented 
in  Figures  27-28  where  comparison  of  the  FEA  and  experimental  results  are  presented. 
The  cohesive  model  parameters  were  chosen  so  that  that  good  agreement  was  obtained 
for  Sample  2.  Thus,  comparing  the  results  of  the  remaining  samples,  the  rate-independent 
FEA  results  don’t  agree  as  suspected,  indicating  a  rate-dependent  model  is  required  to 
adequately  model  this  material  subjected  to  high  loading  rates. 


Figure  25.  Measured  Crack  Opening  at  Specimen  End  versus  time  (Samples  1, 2,  and  4). 
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Figure  26.  Measured  Crack  Opening  at  Specimen  End  versus  time  (Samples  5,  6,  and  9). 
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Table  3.  Prescribed  Velocities  Used  for  FEA  of  RHAMM  DCB  Experiments 


Sample 

Time  (ms) 

Velocity 

(mm/ms) 

0-3.1 

3.538 

1 

3. 1-7.1 

6.968 

7.1- 

4.323 

2 

0- 

1.863 

0-2 

0.619 

4 

24.5 

2.367 

4.5- 

4.864 

0-2 

3.941 

5 

2-3.5 

8.386 

4.5- 

14.484 

0-2 

2.147 

6 

2- 

6.813 

0-4 

1.885 

9 

4- 

6.590 
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Table  4.  Cohesive  element  parameters  for  RHAMM  Experiments 


NCT-350-TR-50 

Parameters 

Rate-Independent 

Rate-Dependent 

uSt  pm 

0.01 

0.01 

Her,  pm 

0.344 

2.0 

Herd,  pm 

30 

2.5 

,  GPa 

0.0200 

0.0200 

X  ,  GPa 

max 

0.0275 

0.0275 

v,  pm/ms 

0.0 

0.0 

ij,  pm/ms 

0.0 

220 

Time  (ms) 


Figure  27.  Comparison  of  DCB  Rate-Independent  FEA  to  experimental  (samples  1, 2,  and  4):  crack  length 
versus  time. 
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Figure  28.  Comparison  of  DCB  Rate-Independent  FEA  (samples  5,  6,  and  9)  to  experimental:  crack  length 
versus  time. 


The  rate-dependent  cohesive  model  (rj^Q)  was  then  used  with  a  single  set  of  parameters 
presented  in  Table  3.  A  comparison  of  the  rate-dependent  FEA  results  to  the  experimental 
results  are  presented  in  Figures  29-30.  Although  the  agreement  was  not  perfect,  much 
better  agreement  was  found  with  the  rate-dependent  model  compared  to  the  rate- 
independent  model. 
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Figure  29.  Comparison  of  DCB  Rate-Dependent  FEA  to  experimental  (samples  1 , 2,  and  4):  crack  length 
versus  time. 
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Figure  30.  Comparison  of  DCB  Rate-Dependent  FEA  (samples  5,  6,  and  9)  to  experimental:  crack  length 
versus  time. 
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5  Conclusions 

The  unique  rate-dependent  cohesive  zone  model  presented  in  this  paper  was  implemented 
using  the  user  defined  cohesive  material  model  of  LS-DYNA  and  was  used  to  model 
published  experimental  mode  I  testing  of  epoxy  and  PEEK  /  carbon-fiber  composites  at 
varied  loading  rates.  With  a  rate-independent  model,  agreement  of  the  published 
experimental  data  and  the  FEA  at  varied  loading  rates  was  obtained  only  by  increasing 
the  critical  strain  energy  release  rate.  The  rate-dependent  model,  with  a  single  set  of 
constants  for  each  composite  material,  good  agreement  between  the  FEA  and  previously 
published  experimental  results  was  found  over  the  loading  rates  considered.  “Stick-slip” 
crack  growth  behavior  observed  experimentally  at  specific  rates  was  also  captured. 

The  rate-dependent  model  developed  during  this  research  is  unique  as  it  contains  a  rate- 
dependent  term  for  the  onset  of  the  opening  ( vii )  of  the  cohesive  zone  combined  with  a 
maximum  traction  rate-dependent  term  (utj).  The  model  is  also  unique  in  that  a  cohesive 
zone  starting  or  yielding  displacement  ( us )  may  be  defined  to  greatly  increase  numerical 
stability. 

To  the  author’s  knowledge  this  is  the  first  time  a  rate-dependent  cohesive  model  has  been 
implemented  into  a  commercial  code  for  both  shells  and  bricks  and  has  been  validated 
with  experimental  data  at  varied  loading  rates.  Additionally  with  the  addition  of  the  yield 
displacement  ( us ),  this  method  was  found  to  be  computationally  robust. 

A  simple  inexpensive  testing  apparatus  was  developed  for  high  rate  DCB  testing.  With 
this  apparatus,  samples  may  be  subjected  to  varied  0.6  to  14.5  mm/s.  FEA  of  experiments 
performed  using  this  apparatus  showed  that  the  rate-dependent  cohesive  model  is 
necessary  to  adequately  model  the  composite  material  tested  (NCT-350-TR-5).  Thus, 
using  this  apparatus  the  material  parameters  needed  for  the  rate-dependent  cohesive  zone 
model  may  be  determined  with  a  small  number  of  experiments  (<10). 

The  rate-dependent  cohesive  zone  model  and  experimental  apparatus  developed  in  this 
project  have  been  shown  to  be  viable  for  use  of  accurately  predicting  the  survivability  of 
multi-sparred,  fuel-filled  composite  wing  structures. 
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6  Recommendations 


For  extension  of  this  work,  it  is  recommended  to  expand  the  mode  I  rate-dependent 
cohesive  element  validated  in  this  work  to  include  mode  II  and  mixed  mode  rate- 
dependent  behavior.  This  model  would  be  validated  for  pure  mode  II  using  simple  double 
cantilever  beam  specimens  and  for  mixed  mode  (mode  I  and  II)  using  real  composite 
wing  joints.  Numerical  results  should  be  matched  with  experimental  results  with  varied 
rates  consistent  with  hydrodynamic  ram.  Multiple  joint  geometries  and  materials  should 
be  considered  for  this  validation. 
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A3.  DCB  High  Rate  Experimental  Apparatus  Drawings 


Figure  A3-3.  Assembly  Drawing  of  DCB  High  Rate  Experimental  Apparatus:  Spring  Mechanism. 
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Table  A3-1:  List  of  Parts 


part 

material 

min  dimension 

bearings 

plates 

AL 

.25x12x12 

subframe 

al 

12x3x4 

framepost 

al 

13xd,75 

bolts 

1/4-28x1  in 

bolts 

1/2-20x1 

springhousing 

flange  /  pipe 

springblock 

al 

2x3x4 

pinblock 

al 

3x4x6 

arms 

steel 

5.5x5x.0625 

slider 

al 

hollow  ball 

bolts 

12-28  x.5in 

bolts 

12-28  xlin 

pullpin 

al 

d. 25x1 .75 

nuts 

size  1-12 

Qty 

3 

3 

2 

8 

18 

20 

1 

1 

1 

1 

1 

10 

3 

1 
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Figure  A3-4.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Base  Plate. 


39 


(5.5000 


Figure  A3-5.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Large  Washer. 
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Figure  A3-6.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Plug. 
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Figure  A3-7.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Framepost. 
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Figure  A3-8.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Turnkey. 
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Figure  A3-9.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Sliding  Pin  Block. 


Figure  A3- 10.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Slider  Plate. 
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Figure  A3-1 1 .  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Pull  Pin. 
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Figure  A3- 12.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  PushRod  Mount. 
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Figure  A3- 14.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Lever  Arm. 
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Figure  A3- 15.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Spring  Block. 
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Figure  A3- 16.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Upright 


51 


Figure  A3- 17.  Drawing  for  DCB  High  Rate  Experimental  Apparatus:  Top. 


A4.  LS-DYNA  User  Cohesive  Model  Code 


*  MAT_U S  ER_DE  F I NE  D_MATER I AL_MODEL  S 
$$ 


$#  mid 

ro 

mt 

lmc 

nhv 

iortho 

ibulk 

1 

1 .5400E-6 

42 

15 

6 

$#  ivect 

ifail 

itherm 

ihyper 

ieos 

1 

1 

0 

-1 

$#  ROFLG 

sfrac 

INTFAIL 

uls 

ulint 

ulend 

syield 

smax 

0.000 

1.000000 

1.00E-5 

3 . 4400E-4 

0.03000 

0.025 

0.0275 

10.0 

$#  etacs 

nubdc 

kp 

ckm2 

u2int 

u2end 

fpen 

0.000 

0.0 

1.0 

80.0  3. 

.  4400E-4 

0.020 

0.01 

subroutine  umat42c (idpart, params, 1ft, lit, fTraction, jump__u, dxdt, 
&  aux, ek, ifail, dtlsiz, crv) 

include  'nlqparm' 
c 

c  Rate-Dependent  Vectorized  Cohesive  Material  Uer  Model 

c  Novemeber  15,  2005 

c 

c  Brian  D.  Choules,  Ph.D 

c  RHAMM  Technologies 

c  332  Skyland  Drive 

c  Bellbrook  OH  45305-8717 

c 

c  Inputs 

c  Mode  I  parameters 

c  uls  =  displacement  for  syield 

c  ulint  =  displacement  at  start  of  damage,  ucr 

c  ulend  =  displacement  at  end  of  damage,  uend 

c  syield  =  yield  stress  at  u=uls 

c  smax  =  maximum  stress  at  ucr  and  dudt=0 

c  sfrac  =  fracture  Stress 

c  etacs  =  viscosity  parameter  for  smax 

c  nubdc  =  viscosity  parameter  for  syield 

c  Kp  =  gain  parameter  for  keeping  nodes  together  when  not 

yielded 

c  Mode  II  and  III  parameters 
c  ckm2  =  stiffness  mode  II  and  III 

c  u2int  =  mode  II  and  III  critical  displacment 

c  u2end  =  mode  II  and  III  end  displacement 

c  Other  Input 

c  fpenalty  =  constant  for  resisting  penetration 

c 

c***  variables 

c  idpart  =  part  ID 

c  params  =  material  constants 

c  1ft, lit  =  start  and  end  of  block 

c  fTraction  =  components  of  the  cohesive  force 

c  jump_u  =  components  of  the  displacement 
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c  dxdt  =  components  of  the  velocity 

c  aux  =  history  storage 

c  aux(i,l)  =storage  for  yielding  condition  (yieldd) 

c  aux(i,2)  =storage  for  last  time  step  jump  (dinl) 

c  aux(i,3)  =storage  for  critical  crack  opening  rate  at  ucr 

(dudtcrtemp) 

c  aux(i,4)  =storage  for  Damage  Variable  (daml) 

c  aux(i,5)  =storage  for  last  stress  for  keeping  nodes 

together  when  us=0  (signl) 

c  aux(i,6)  =storage  for  critical  crack  opening  rate  at  us 

(dudtstemp) 

c  ek  =  max.  stiffness/area  for  time  step  calculation 

c  ifail  =. false,  not  failed 

c  =.true.  failed 

c  dtlsiz  =  time  step  size 

c  crv  =curve  array  -  not  used 

c  dil  =  normal  displacement  jump 

c  di2  =  mode  II  displacement  jump 

c  di3  =  mode  III  displacement  jump 

c 

c  daml  =  mode  I  damage  variable 

c  di23  =  mode  II  and  III  combined  damage  variable 

c  dudt  =  mode  I  crack  opening  rate 

c 

£  ********** 

c 

c  jump_u,  dxdt,  and  fTraction  are  in  the  local  coordinate  system: 
c  components  1  and  2  are  in  the  plane  of  the  cohesive  surface 
c  component  3  is  normal  to  the  plane 
c 

c  the  declaration  below  is  processed  by  the  C  preprocessor  and 

c  is  real*4  or  real*8  depending  on  whether  LS-DYNA  is  single  or 

double 

c  precision 

real  L, jump_u, cklb, ck2b, uls, ulint, u2int , ulend,  u2end, 

&  ul delta, u2delta, ulie, u2ie, daml , damli, dam23i, dam23 
&  ckm2, dudt, dudtcrtemp, dudtcrl , dudtstemp, dudts 
&  etacs, alphabdc, syield, smax, sfrac, nubd 
&  fK, fac, fpenalty 

&  Kp,  dil , di2, di3, di23, dinl , din2, din3, signl , one, zero 
logical  ifail 

dimension  params (* ) , fTraction (nlq, *) , jump_u (nlq, * ) , 

&  dxdt (nlq, *) , aux (nlq,*),ek(*) , ifail (*) ,dtlsiz(*) , 

&  crv (101,2,*) 

c  connector  stiffnesses/area  (3  direction  normal  to  delam  plane) 

c 

c  critical  separation  displacement  to  fail  connector 

one=l .  0 
zero=0 . 0 

c  Mode  I  parameters 

uls  =  params (3)  !  displacement  for  syield 

ulint  =  params (4)  !  displacement  at  start  of  damage,  ucr 

ulend  =  params (5)  !  displacement  at  end  of  damage,  uend 

syield  =  params (6)  !  yield  stress  at  u=uls 
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smax  =  params(7)  !  maximum  stress  at  ucr  and  dudt=0 
sfrac  =  params(8)  !  fracture  Stress 
etacs  =  params(9)  !  viscosity  parameter  for  smax 
nubdc  =  params(lO)  !  viscosity  parameter  for  syield 
Kp  =  params(ll)  !  gain  parameter  for  keeping  nodes  together  when 
not  yielded 

c  Mode  II  and  III  parameters 

ckm2  =  params(12)  !  stiffness  mode  II  and  III 
u2int  =  params(13)  !  mode  II  and  III  critical  displacement 
u2end  =  params(14)  !  mode  II  and  III  end  displacement 
c  Other 

fpenalty  =  params(15)  !  penetration  constant 
c 

fK=fpenalty*smax/ulint 
c  intermediate  material  constants 
uldelta  =  ulint  -  ulend 
u2delta  =  u2int  -  u2end 
ulie  =  ulint*ulend 
u2ie  =  u2int*u2end 
c  Main  do  loop 

do  i=lft,llt 

di2  =  jump_u(i,l)  !  mode  II 
di3  =  jump_u(i,2)  !  mode  III 
dil  =  jump_u(i,3)  !  normal  displacement 
c 

di23=sqrt (di2*di2+di3*di3)  !  combined  mode  II  and  III 
displacement 

c  Extract  History  Variables 

yieldd=aux (i, 1)  !  element  yielded  when  equal  one  -  not 
yeilded  when  qual  to  zero 

if (uls.eq. zero. and. syield. eq. zero)  then 
yieldd=one 
aux (i, l)=one 
endif 
c 

dinl=  aux(i,2)  !  last  time  step  jump 
aux(i,2)=dil  !  set  for  next  time  step 
dudtcrl=aux (i, 3)  Hast  time  step  dudtcrtemp 
dudts=aux (i, 6)  Hast  time  step  dudtcrtemp 

signl=aux (i, 5)  !  last  stress  for  keeping  nodes  together  prior 

to  yield 
c 

c  Set  damage  variables 
c 

c  Mode  I  damage  variable 

if (dil . It .ulint)  then 
daml  =  one 
aux(i,4)=  one 

else  if (dil .ge. ulint)  then 

damli  =  (dil*ulint-ulie) / (dil*uldelta)  !  calculate 
intermediate  damage  variable 

daml  =  max (damli, zero)  !if  negative  set  to  zero 

daml  =  min (aux (i, 4) , daml)  !  store  minimum  positive  value 

aux(i,4)  =  daml  (store  it 

if (daml . le . zero)  ifail (i)=.true. 

1  endif 

c  Mode  II  and  III  damage  variable 
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if (di23.gt.u2int) then 

dam23i=  (di23*u2int-u2ie) / (di23*u2delta)  !  calculate 
intermediate  damage  variable  mode  II  and  III 

dam23=max (dam23i, zero)  !  if  negative  set  to  zero 
if (dam23 .le.zero)  ifail  (i)  =  . true . 
else 

dam23=one 

endif 

c 

C  CALCULATE  STRESSES 

c 

c  Rate  dependent  term  -  dudt 
c 

c  Store  cricital  dudt  value  -  Last  value  that  gets  used  at  ulint 
if(dil  .le.  uls)  then 

if ( (dil-dinl )  .gt.  zero  .and.  dtlsiz(i)  .gt.  zero  )  then 
dudt= ( (dil-dinl ) /dtlsiz (i) ) 
else 

dudt=zero 

endif 

dudtstemp=  max (dudts, dudt)  !new  or  old  value 
aux(i,6)  =  dudtstemp  (store  it  critical  dudt 
else  if(dil  .le.  ulint)  then 

if ( (dil-dinl)  .gt.  zero  .and.  dtlsiz (i)  .gt.  zero  )  then 
dudt= ( (dil-dinl) /dtlsiz (i) ) 
else 

dudt=zero 

endif 

dudtcrtemp=  max (dudtcrl, dudt)  (new  or  old  value 
aux(i,3)  =  dudtcrtemp  (store  it  critical  dudt 

else 

dudtcrtemp=dudtcrl 

dudtstemp=dudts 

endif 

c  Assign  Tractions 

if (uls . eq. zero) then 

if (yieldd.eq.zero)  then 
if (dil . le . zero)  then 

fTraction (i, 3) =  dil*syield/ulint+fK*dil 
else 

fTraction (i, 3) =  dil* (smax/ulint) *Kp+signl 
if  ( fTraction (i, 3) . ge . syield)  then 
aux (i, 1) =one 
yieldd=one 
endif 
endif 

else  !  already  yielded 

if  (ifail (i) .eq. true)  then 
fTraction (i, 3) =0 . 0 
else  if (dil . le . zero)  then 

fTraction (i, 3) =  dil*syield/ulint+f K*dil 
else  if (dil. le. ulint)  then 

fTraction (i, 3) =syield+nubdc*dudtstemp+ (smax-syield  + 
&  dudtcrtemp*etacs-dudtstemp*nubdc) * ( (dil-uls) / (ulint- 

uls)  ) 

c  check  for  negative  tractions  which  can  occur  with 

etacs  negative 
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if (fTraction (i, 3) .It . zero)  then 
f Tract ion (i, 3) =0 . 0 
ifail (i) =.true. 
endif 
else 

fTraction (i, 3) =  daml*dil* 

&  (smax+etacs*dudtcrtemp) /ulint 

c  check  for  negative  tractions  which  can  occur  with 

etacs  negative 

if (fTraction (i, 3) . It . zero)  then 
fTraction (i, 3) =0.0 
ifail  (i)  = . true . 
endif 
endif 
endif 

else  !  uls  not  zero 

if  (ifail (i) .eq. true)  then 
fTraction  (i, 3) =0 . 0 
else  if (dil . le . zero)  then 

fTraction (i, 3) =  dil*syield/uls+fK*dil 
else  if (dil .le. uls)  then 

fTraction (i, 3) =  (syield+nubdc*dudtstemp) *dil/uls 
c  check  for  negative  tractions  which  can  occur  with  etacs 

negative 

if (fTraction  (i, 3) . It . zero)  then 
fTraction (i, 3) =0. 0 
ifail  (i)  =  .true. 
endif 

else  if (dil .le. ulint)  then 

fTraction  (i, 3) =syield+nubdc*dudtstemp+ (smax-syield  + 

&  dudtcrtemp*etacs-dudtstemp*nubdc) * ( (dil-uls) / (ulint- 

uls)  ) 

c  check  for  negative  tractions  which  can  occur  with  etacs 

negative 

if (fTraction (i, 3) . It . zero)  then 
fTraction (i, 3) =0 . 0 
ifail (i)=.true. 
endif 
else 

fTraction (i, 3) =  daml*dil* 

&  (smax+etacs*dudtcrtemp) /ulint 

if (fTraction (i, 3) . It . zero)  then 
fTraction  (i, 3) =0 . 0 
ifail  (i)  =  .true. 
endif 
endif  " 

c  YIELD  criteria  when  syield  >0 

if (dil .ge .uls . and. yieldd.eq. zero .and. syield .gt .zero)  then 
aux(i,l)=one  !  store  yieldd  in  history  variable 

endif 
endif 
c 

c  Mode  2  and  Mode  3  no  rate  dependence 

fTraction (i, 1) =ckm2*di2*dam23 
fTraction (i, 2) =ckm2*di3*dam23 
c 

c  approximate  stiffness  for  timestep 
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if  (dil  .It.  0)  then 

ek  (i) =f Tract ion (i, 3) /dil  +  fK 
if (ek(i) .gt. (smax*1000000. /ulint) )then 
write  (13,  1050) 
write  (*,  1050) 
ek (i) =smax* 1000000 . /ulint 
endif 
else 

ek (i) =f Tract ion (i, 3) /dil 
if (ek (i) .gt . ( smax* 1000000 . /ulint)  )  then 
write (13, 1060) 
write (*, 1060) 
ek (i) =smax*1000000 . /ulint 
endif 
endif 
c 

c  FAILURE  CRITERIA 

if (f Tract ion (i, 3) .gt . sf rac .and. if ail (i) .eq.  .false. 

&  . and. sf rac .ne . zero)  then 

if ail (i) = . true .  !  for  failure  criteria 
aux(i,l)=one  !  store  yield  in  history  variable 

c  diognostic  info 

write  (13, 1010)  i, fTraction (i, 3) 
write (*, 1010)  i, fTraction (i, 3) 
endif 

enddo  lend  of  the  i  loop 

c 

return 

1010  format (5x, ' delamination  element  has  stress  greater  than  sfrac'i8, 
&  lpel2 . 4 ) 

1050  format (5x, 'ek  limit  reached  on  penetration') 

1060  format (5x, ' ek  limit  reached') 
end 
c 
c 
c 
c 
c 
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A5.  Critical  Portions  LS-DYNA  Input  File  for  DCB  Model 

* KEYWORD 
*TITLE 

EPOXY  CARBON  FIBER 
*  CONTROL_PARALLE  L 
4 

* CONTROL  BULK  VISCOSITY 


$#  ql 

q2 

type 

1.500000 

0.060000 

‘CONTROL  CONTACT 

$#  slsfac 

rwpnal 

islchk 

enmass 

0.100000 

0.000 

1 

$#  usrstr 

usrf rc 

nsbcs 

tiedpr j 

0 

0 

10 

$#  sfric 

df  ric 

edc 

0.000 

0.000 

0.000 

$#  ignore 

f rceng 

skiprwg 

0 

0 

0 

‘CONTROL  COUPLING 

$#  unleng 

untime 

unforc 

subcyl 

1.000000 

1.000000 

1.000000 

1 

*CONTROL_CPU 
$#  cputim 


0.000 

‘CONTROL  DYNAMIC  RELAXATION 


$#  nrcyck 

drtol 

drfctr 

idrf lg 

250 

0.001000 

0.9950001 

‘CONTROL  ENERGY 

$#  hgen 

rwen 

slnten 

1 

2 

1 

‘CONTROL  HOURGLASS 

$#  ihq 

qh 

1 

0.100000 

‘CONTROL  OUTPUT 

$#  npopt 

if lush 

neecho 

nrefup 

0 

3 

0 

$#  iprtf 

0 

‘CONTROL  SHELL 

$#  wrpang 
proj 

esort 

irnxx 

20.000000 

0 

-1 

$#  rotascl 

intgrd 

lamsht 

1.000000 

0 

0 

‘CONTROL  TERMINATION 

$#  endtim 

endcyc 

dtmin 

2.00000 

0 

0.000 

‘CONTROL  TIMESTEP 

shlthk 

penopt 

thkchg 

1 

1 

1 

interm 

xpene 

ssthk 

0 

4.000000 

vfc 

th 

th  sf 

0.000 

0.000 

0.000 

outseg 

spotstp 

spotdel 

0 

0 

0 

timidl 

flipx 

flipy 

0.000 

0 

0 

drterm 

tssfdr 

irelal 

OOOOE+30 

0.900000 

0 

rylen 

1 

iaccop 

opif  s 

ipnint 

0 

0.000 

istupd 

theory 

bwc 

0 

cstyp6 

1 

2 

tshell 

2 

nfaill 

endeng 

0.000 

endmas 

0.000 

orien 

1 

ecdt 

pen_sf 

0.000 

f  lipz 
0 

edttl 

0.040000 


ikedit 

miter 

1 

nf ai!4 
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$#  dtinit 
msist 

tssfac 

isdo 

tslimt 

dt2ms 

lctm 

erode 

0.000 

0.670000 

0 

0.000 

0.000 

$#  dt2msf 

dt2mslc 

0.000 

*  DAT ABAS  E_B I NARY_D 3  PLOT 

$#  dt  lcdt  beam  npltc 

0.02 

$#  ioopt 

0 

*BOUNDARY_PRESCRIBE  D_MO  T I ON_S E T_ I D 

$#  cid 

heading 

lUpper  Arm  Velocity 

$#  nsid  dof  vad  lcid  sf  vid  death 

birth 

2401  20.50000  2  100.0000 

0.000 

$* BOUNDARY  JSPC_SET_ID 

$#  cid 

heading 

$  lB'ottom  Arm  Fixed 

$ 

$  LBC  set  :  Fixed  Beam 
$ 


$# 

dof  rz 

nsid 

cid 

dofx 

dofy 

dof  z 

dof  rx 

dofry 

$ 

1 

0 

0 

1 

0 

0 

0 

yj 

*SET  NODE  LIST 

TITLE 

Fixed 

BEam  2 

$# 

sid 

dal 

da2 

da3 

da  4 

1 

0.000 

0.000 

0.000 

0.000 

$# 

nid8 

nidi 

nid2 

nid3 

nid4 

nid5 

nid6 

nid7 

11623 

11624 

11625 

11626 

11627 

11628 

11629 

11630 

$  nodes  left 

out  for  brevity 

$$$ 

$* CONTACT  AUTOMATIC  GENERAL  ID 

$# 

title 

cid 

$ 

lcontact  DCB 

$# 

ssid 

msid 

sstyp 

mstyp 

sboxid 

mboxid 

spr 

mpr 

$ 

0 

$# 

0 

0 

0 

0 

0 

0 

0 

f  s 

fd 

dc 

VC 

vdc 

penchk 

bt 

dt 

$ 

0.200 

0.000 

0.000 

0.000 

0.000 

0 

0.0001 

. 0000E+20 

$# 

sf  s 

sfm 

sst 

mst 

sf  st 

sfmt 

fsf 

vsf 

$  1.000000 
1.000000 

1. 

000000 

0.000 

0.000 

1.000000 

1.000000 

1.000000 

*PART 

$#  title 

60 


CZPROP 

$#  pid  secid  mid  eosid  hgid  grav  adpopt 

tmid 

1  1  1 
*  SECT I ON_SOL I D_T I TLE 
Section  Cohesive  Zone 


$#  secid 

elform 

aet 

1 

20 

1 

$$ 

*MAT  USER  DEFINED  MATERIAL  MODELS 

$$ 

$#  mid 

ig 

ro 

mt 

lmc 

nhv 

iortho 

ibulk 

1  1 

. 5400E-6 

42 

15 

6 

$#  ivect 

ifail 

itherm 

ihyper 

ieos 

1 

1 

0 

-1 

$#  ROFLG 

sf  rac 

INTFAIL 

uls 

ulint 

ulend 

syield 

smax 

0.000 

1.000000 

1.00E-5 

3 . 4400E-4 

0.00700 

0.02 

0.0275 

10.0 

$#  etacs 

nubdc 

kp 

ckm2 

u2int 

u2end 

fpen 

0.05 

0.0 

1.0 

80.0  3. 

.  4400E-4 

0.020 

0.01 

*PART 
$#  title 
PROPSHELL 
$#  pid 

tmid 

secid 

mid 

eosid 

hgid 

grav 

adpopt 

2 

2 

2 

* SECT I ON  SHELL  TITLE 

Section  DCB 

Beam 

$#  secid 

setyp 

elform 

shrf 

nip 

propt 

qr/irid 

icomp 

2 

1 

16 

0.830 

5 

1 

0.000 

0 

1 

$#  tl 

t2 

t3 

t4 

nloc 

marea 

1.500000 

1.500000 

1.500000 

1.500000 

0 

0.000 

*MAT  ELASTIC 

TITLE 

Shell  Composite 

$#  mid 

ro 

e 

Pr 

da 

db 

not  used 

2  1 

.  5660E-6 

120.00000 

0.270000 

0.000 

0.000 

* INITIAL_VELOCITY 
$ 

$  LBC  set  :  initVelocity 
$ 


$# 

nsid 

3 

nsidex 

boxid 

irigid 

$# 

vx 

vy 

vz 

vxr 

vyr 

vzr 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

*DEFINE_CURVE_TITLE 
Upper  Arm  Velocity  Curve 
$ 

$  This  is  the  unit  load  curve 
$ 


$# 

lcid 

sidr 

sfa  sfo 

of  fa 

of  fo 

dattyp 

1 

0 

0.000  0.000 

0.000 

0.000 

$# 

al 

ol 

0.000 

1.00000000 
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1.00000000 


100.00000000 
*DEFINE__VECTOR 
$ 

$  LBC  set  :  Velocity  UpperBeam 


$ 

$#  vid 

xt 

yt 

1 

0.000 

0.000 

*DEFINE_VECTOR 
$#  vid 

xt 

yt 

2 

0.000 

0.000 

*DEFINE_VECTOR 
$#  vid 

xt 

yt 

3 

0.000 

0.000 

*SET  NODE  LIST 

TITLE 

Velocity  Contorl  NODES 

$#  sid 

dal 

da2 

2 

0.000 

0.000 

$#  nidi 

nid8 

nid2 

nid3 

12403 

12404 

12405 

12410 

$  nodes  left  out  for  brevity 

*SET  NODE  LIST 

TITLE 

ALL  NODES  for  Initial 

Velocity 

$#  sid 

dal 

da2 

3 

0.000 

0.000 

$#  nidi 

nid8 

nid2 

nid3 

1 

2 

3 

8 

$  nodes  left  out  for  brevity 
*DAMPING_PART_STIFFNESS 
2  0.01 
*ELEMENT_SOLID 

$#  eid  pid  nl  n2 

n7  n8 

111  443 

665  223 

$  nodes  left  out  for  brevity 
*ELEMENT_SHELL 

$#  eid  pid  nl  n2 

5701  2  11493  11494 

$  element  definitions  left  out  for 
*NODE 

$#  nid  x 

rc 

1  30 

$  nodes  left  out  for  brevity 
*END 


zt 

0.000 

xh 

1.000000 

yh 

0.000 

zh 

0.000 

zt 

0.000 

xh 

0.000 

yh 

1.000000 

zh 

0.000 

zt 

0.000 

xh 

0.000 

yh 

0.000 

zh 

1.000000 

da  3 
0.000 
nid4 

da  4 
0.000 
nid5 

nid6 

nid7 

12406 

12407 

12408 

12409 

da3 

0.000 

nid4 

da  4 
0.000 
nid5 

nid6 

nid7 

4 

5 

6 

7 

n3 

n4 

n5 

n6 

444 

2 

222 

664 

n3 

11520 

brevity 

n4 

11519 

y 

z 

tc 

-0.75 

0 
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A6.Crack  Growth  Data  Extraction  Code 


/*  This  program  takes  curve  file  from  Ispost  with  nodes  on  both  sides 
of  an  interface  and  outputs  a  vs  t  */ 

/*  The  crack  is  assumed  to  be  in  the  XZ  plane  with  the  crack  length 
growth  direction  in  the  X  direction  */ 

/*  the  crack  length  (a)  is  defined  as:  distance  from  edge  to  crack  top 
minus  distance  from  edge  to  loading  pin  location  (edist)  */ 

/*  A  file  with  node  xyz  coordinates  must  be  specified  in  following 
format 
♦NODE 

1  30.00000000  -0.75000000  0.000 


endf ile 

*/ 

#include  <math.h> 

♦include  <stdio.h> 

♦include  <stdlib.h> 

♦include  <string.h> 

♦define  NMAX  30000  /*  Intial  guess  at  Maximum  number  of  NODES  */ 

♦define  ELEMMAX  25000  /*  Intial  guess  at  Maximum  number  of  ELEMENTS  */ 

♦define  LINESMAX  3000000  /*  MAX  number  of  lines  to  scan  */ 

♦define  TOL  0.0011 

♦define  NPOPEN  1  /*  Number  of  consecutive  nodal  locations  where  the 
crack  opening  must  be  greater  than  Umax  */ 

int  main  () 

{ 

int  numxypairs,  numnodepairs, garbage; 
int  nodenumtemp; 
int  i ,  j  ,  k ; 

int  numnodes;  /*  total  number  of  nodes  */ 

int  *nodenum;  /*  matrix  of  node  and  element  numbers  (IDS)*/ 
double  *plottime; 

double  **time,  **ydisp;  /*time  and  displacement  matrices  */ 
double  *acrack; 
double  dydisp; 

double  aO;  /*  Initial  crack  length  from  load  point  */ 

double  edist;  /*  Distance  from  End  of  sample  to  load  location  */ 

int  oldNM=NMAX; 

int  newNM=NMAX; 

int  *nodesearch; 

int  **nodematch; 

char  str[20]; 
char  filename [20 ) ; 
char  filename2 [20] ; 

char  str2 [200] ; 


63 


FILE  *fp,  *fn; 
double  *x,*y, *z; 

double  xtemp, ytemp, ztemp; 

double  umax;  /*  maximum  Distance  between  cohesive  element  nodes*/ 
void  qsortl (double  *a,  double  *b,int  n) ; 


/*******Allocate  Matrices  ****************/ 

/*  node  matrices  */ 

nodenum=malloc (NMAX*sizeof (int) ) ; 

/*  Read  input  file  */ 

printf ("Enter  complete  filename\n") ; 

gets (str) ; 

sscanf (str, "%s" , filename) ; 
fp  =  fopen (filename, "r") ; 
if ( fp  ==  NULL) 

{ 

printf ("file  not  available\n") ; 
exit  (1) ; 

} 

/*  Enter  in  umax  */ 

printf  ( "Enter  umax (mm)  \n")  ; 

gets  (str) ; 

sscanf (str, "%lf", Sumax) ; 

/*  Enter  in  aO  */ 

printf  ( "Enter  ao(mm)  \n"),- 

gets (str) ; 

sscanf (str, "%lf", &a0) ; 

/*  Enter  in  e  */ 

printf  ("Enter  edist  (mm)  \n"),- 

gets (str) ; 

sscanf (str, "%lf", Sedist) ; 

/*  looking  for  Node  no.  Keyword  */ 
for (i=0; i<=LINESMAX; i++) { 

fgets  (str2, sizeof (str2) , fp) ; 
if (strncmp (str2, "Node  no . 8) ==0) { 
break; 

} else  if (strncmp (str2, "endcurve", 8) ==0) ( 
printf ("endcurve  found\n") ; 
exit  (1) ; 

) 

} 


fgets (str2, sizeof (str2) , fp) ; 

sscanf (str2 , "%d  %5s  %d", &nodenum[0] , str,  Snumxypairs) ; 
/*  Allocate  time  and  ydisp  */ 
time=malloc (NMAX*sizeof (double) ) ; 
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ydisp=malloc (NMAX*sizeof (double) ) ; 
for (i=0; i<NMAX; i++) 

( 

time [i]=malloc (numxypairs*sizeof (double) ) ; 
ydisp[i]=malloc (numxypairs*sizeof (double) )  ; 

) 

i=0; 

/*  Reading  in  time  disp  for  each  node*/ 
while (strpbrk (str2, "#pts=") !=NULL) { 
for (j=0; j<numxypairs; j++) { 

fgets (str2, sizeof (str2) , fp) ; 

sscanf (str2, "%lf  %lf  ", stime [i] [ j ] , Sydisp [i] [ j ] ) ; 

} 

i++  ; 

/*  Reallocate  Matrices  for  time  and  x  y  data  */ 
if (i>=newNM-l) ( 

printf ("im  am  reallocing  node  matrices  "); 

oldNM=newNM; 

newNM= (int) (oldNM*1.2); 

nodenum  =  realloc (nodenum, newNM*sizeof (int  *)); 
if  (nodenum==NULL)  (printf("not  enough  memoryXn"); 

exit  (1)  ;  } 

time  =  realloc (time,  newNM*sizeof (double  *)); 
if  (time==NULL)  (printf ("not  enough  memoryXn"); 

exit  (1)  ;  ) 

ydisp  =  realloc (ydisp, newNM*sizeof (double  *)); 
if  (ydisp==NULL)  (printf ("not  enough  memoryXn"); 

exit  (1) ;  } 

for  (i=oldNM; i<newNM; i++) 

{ 

time[i]  =  malloc (numxypairs*sizeof (double) ) ; 
if  ( t ime [ i ) ==NULL )  ( 

printf  ("not  enough  memoryXn");  exit(l);} 
ydisp[i]  =  malloc (numxypairs*sizeof (double) ) ; 
if  (ydisp[i]==NULL)  ( 

printf ("not  enough  memoryXn");  exit(l);} 

} 

} 

fgets (str2, sizeof (str2) , fp) ;  /*  get  endcurve  out  of  the  way 

*/ 

fgets (str2, sizeof (str2) , fp) ;  /*  get  Node  #  and  #  points  */ 
sscanf (str2, "%d  %6s  %d", Snodenum [i] , str,  Sgarbage) ; 

} 

numnodes=i; 
f close ( fp) ; 

x=malloc (newNM*sizeof (double)  )  ; 
y=malloc (newNM*sizeof (double)  )  ; 
z=malloc (newNM*sizeof (double) ) ; 

! *  **************************************************** j 

/*  Read  node  coordinate  input  file  */ 
printf ("Enter  node  coordinate  filenameXn") ; 
gets (str)  ; 
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sscanf  (str, "%s", filename2)  ; 


fp  =  fopen {filename2, "r") ; 
if (fp  ==  NULL) 

{ 

printf ("file  not  available\n") ; 
exit  (1)  ; 

} 

/*  READ  IN  NODE  COORDINATES  */ 

/*  looking  for  *NODE  Keyword  */ 
for  (i=0; i<=LINESMAX; i++) < 

fgets (str2, sizeof  (str2) , fp) ; 
if (strncmp(str2, "*NODE", 5) ==0) { 

printf ("*NODE  found  :-}\n"); 
break; 

}else  if (strncmp (str2, "endfile", 7) ==0) { 
printf ("endfile  found\n"); 
exit  (1)  ; 

) 

) 

printf ("numnodes  =  %d\n" , numnodes) ; 

/*  Reading  in  X,Y,Z  coordinates  of  nodes*/ 
i=0; 

fgets (str2, sizeof (str2) , fp) ; 
while (strncmp (str2, "endfile", 7) !=0) { 
sscanf (str2, "%d  %lf  %lf 
%lf " , Snodenumtemp, Sxtemp, sytemp, &ztemp) ; 

for ( j=0; j<numnodes; j++) { 

if (nodenumtemp==nodenum [ j] ) { 
x [ j ]=xtemp; 
y[j]=ytemp; 
z [ j ] =ztemp; 
break; 

} 

} 

fgets  (str2, sizeof (str2) ,  fp)  ; 
i++; 

} 

fclose (fp) ; 

/****★**********************★*******************■***********★****★*****/ 

/***************************************★****★***★***************★*★/ 

/*  Allocate  matrices  for  matching  */ 
nodesearch=malloc (numnodes*sizeof (int) )  ; 

for (i=0; i<numnodes; i++)  { nodesearch [i] =0; ) 

nodematch=malloc ( (int) (numnodes/2) *sizeof (int) ) ; 
for (i=0;i< ( (int)  (numnodes/2) )  ;i+  +  ) 

{ 

nodematch [ i] =malloc (4*sizeof (int) )  ; 

) 

/*  Find  coincident  NODES  */ 
k=0  ; 

for (i=0; i<numnodes; i++) ( 

if  (nodesearch [i]==l) ( 
continue; 


66 


}else{ 

for ( j=i+l ; j<numnodes; j++) { 

if ( (fabs ( (x [ j ] -x [i] ) ) <=TOL) && (fabs ( z [ j ] — 

z  [i] ) <=TOL) ) { 

nodematch [k] [0]=nodenum[i] ; 

nodematch [k] [ 1 ] =nodenum [ j ] ; 

nodematch [k] [2]=i; 

nodematch[k] [ 3 ) = j ; 

nodesearch [ j ] =1 ; 

k++; 

break; 

} 

} 

if (nodesearch [j ] ==0) { 

printf ("Could  not  find  match  for  node 

%d\n",nodenum[i] ) ; 

exit  ( 1 )  ; 

} 

} 

} 

numnodepairs=  (int)  (numnodes/2) ; 

if (k ! =numnodepairs) { 

printf ("Not  everything  matched  !\n"); 
exit  (1 ) ; 

} 

/♦allocate  plottime  and  acrack  */ 
plottime=malloc (numnodepairs*sizeof (double) ) ; 
acrack=malloc (numnodepairs*sizeof (double) ) ; 

/*  Calculated  difference  in  ydisplacements  */ 
for (i=0; i<numnodepairs; i++) { 

for (j=0; j<numxypairs; j++) { 

if ( (y (nodematch [i] [2] ) -y [nodematch (i) [3] ] ) >0 . 0) { 
dydisp=ydisp [nodematch [i]  [2] )  [ j ]  — 
ydisp [nodematch [i]  [3] ]  [  j ]  ; 

}else( 

dydisp=ydisp [nodematch [i]  [3] ]  [j]- 
ydisp [nodematch [ i ]  [2] ]  [  j  ]  ; 

} 

if (dydisp>umax) ( 

plottime [ i ] =time [ i ] [ j ] ; 
if ( (x [nodematch [ i ] [2] ] -edist) <a0) { 
acrack[i]=aO; 

}else( 

acrack [ i] =x [nodematch [i] [2] ] -edist; 


} 

break; 

>else{ 

plottime [ i] =-5 . 0; 

! 


} 

} 

/♦Put  in  Time  order  */ 

qsortl (plottime,  acrack,  numnodepairs) ; 
fn  =  fopen (strcat  (filename, "_plot") , "w”) ; 
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fprintf (fn, "time (ms) \ta (mm) \n") ; 
/*  Put  results  in  _plot  file  */ 
for (i=0; i<numnodepairs; i++) { 
if (plot time [ i] ! =-5 . 0) ( 

fprintf (fn, "%3 . 91f\t 
%5 . 21f \n", plot time [ i] , acrack [i] )  ; 

) 

) 

fclose (fn) ; 
return  (0)  ; 

) 
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Symbols 


a 

— 

mode  I  damage  term  offset  constant 

«23 

= 

mode  I  damage  term  offset  constant 

do 

= 

initial  crack  length  of  DCB  specimen 

b 

= 

mode  11  damage  term  slope  constant 

b2} 

= 

mode  II  damage  term  slope  constant 

Cs 

= 

characteristic  wave  speed  (- ) 

G" 

= 

mode  I  critical  strain  energy  release  rate 

h, 

= 

arm  thickness  of  DCB 

h2 

= 

arm  thickness  of  DCB 

Kn 

= 

cohesive  model  hardening  stiffness 

t 

= 

time 

u 

= 

mode  I  crack  opening  displacement 

un 

= 

mode  II  crack  opening  displacement 

um 

= 

mode  III  crack  opening  displacement 

M23 

= 

root  mean  square  combined  mode  II  and  II  loading 

it 

= 

rate  of  mode  I  crack  opening  displacement 

Ucr 

= 

critical  mode  I  crack  opening  displacement 

M23  cr 

= 

critical  mode  II/II1  crack  opening  displacement 

"cr 

= 

historical  mode  1  rate  of  crack  opening  displacement  when  u=Ucr 

"end 

= 

mode  I  crack  opening  displacement  where  tractions  are  reduced  to  zero 

M23«i  d 

= 

mode  II/III  crack  opening  displacement  where  tractions  are  reduced  to  zero 

Us 

= 

yield  mode  1  crack  opening  displacement 

= 

historical  mode  I  rate  of  crack  opening  displacement  when  u=us 
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1 

p 

V 


p 

T 


''111 


T 


max 


w 


^23 


cohesive  zone  viscosity  for  shift  of  fm2X 
shear  modulus 

cohesive  zone  viscosity  for  shift  of  T 
density 

mode  I  cohesive  traction 
mode  II  cohesive  traction 
mode  III  cohesive  traction 
critical  cohesive  traction 
rate-independent  model  constant 
maximum  cohesive  traction 
mode  I  yield  cohesive  traction 
mode  I  damage  variable  of  cohesive  zone 
damage  variable  of  cohesive  zone 


70 


Acronyms 


CAI  - 

Composites  Affordability  Initiative 

COD- 

Crack  Opening  Displacement 

CMOD  - 

Crack  Mouth  Opening  Displacements 

DCB  - 

Double  Cantilever  Beam 

DFC- 

Air  Force’s  Decoupled  Fuel  Cells 

DLL- 

Design  Limit  Load 

FEA- 

Finite  Element  Analysis 

HEI  - 

High  Explosive  Incendiary 

HRAM- 

Hydrodynamic  Ram 

JLF- 

Joint  Live  Fire 

LFT&E  - 

Fire  Testing  and  Evaluation 

M&S  - 

Modeling  and  Simulation 

71 


